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The study of liquid crystals at equilibrium has led to fundamental insights into the na¬ 
ture of ordered materials, as well as to practical applications such as display technologies. 
Active nematics are a fundamentally different class of liquid crystals, driven away from 
equilibrium by the autonomous motion of their constituent rod-like particles®. This in¬ 
ternally generated activity powers the continuous creation and annihilation of topological 
defects, which leads to complex streaming flows whose chaotic dynamics appear to destroy 
long-range order 5 n . Here, we study these dynamics in experimental and computational 
realizations of active nematics. By tracking thousands of defects over centimetre-scale dis¬ 
tances in microtubule-based active nematics, we identify a non-equilibrium phase charac¬ 
terized by system-spanning orientational order of defects. This emergent order persists 
over hours despite defect lifetimes of only seconds. Similar dynamical structures are ob¬ 
served in coarse-grained simulations, suggesting that defect-ordered phases are a generic 
feature of active nematics. 

Topological defects play important roles in diverse phenomena ranging from high-energy 
physics and cosmology to traditional condensed matter systems 12 -'. For instance, the sponta¬ 
neous unbinding of dislocation pairs mediates the melting of 2D crystals 13 . Despite their usual 
role as centers of disorder, defects can also organize into higher-order equilibrium structures 
with emergent properties, such as liquid crystalline twist-grain-boundary phases and flux-line 
lattices in superconductor^? 15 . Far less is understood about the role of defects in active matter 
systems, which are driven away from equilibrium by the motion of their constituent particles 
ESH 24 . Previous work on active nematics has demonstrated an instability at large wavelengths 5 
which leads to spontaneous defect nucleation and unbinding® 0 . In contrast to the well-studied 
passive defects found in equilibrium matter, defects in active nematics are motile 25 , and are 
continuously generated and annihilated, producing a dynamical defect-riddled phase that is 
inherently non-equilibrium. The observed dynamics are complex and chaotic, and appear to de¬ 
stroy the long-range ordering of the underlying nematic. Here, by tracking thousands of defects 
over long times, we demonstrate that defects self-organize into a higher-order phase with bro¬ 
ken rotational symmetry. The orientational ordering of defects spans macroscopic samples and 
persists for the sample lifetime of many hours, despite the lifetimes of the constituent defects 
being orders of magnitude shorter. 

Our experimental system is comprised of micron-long stabilized microtubules (MTs), 
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streptavidin clusters of biotin-labeled kinesin motors- 6 and the non-adsorbing polymer polyethy¬ 
lene glycol (PEG) (Fig. la). In a bulk suspension, PEG induces formation of MT bundles by the 
depletion mechanism 27 28 . The same interaction also depletes MTs onto a surfactant-stabilized 
oil-water interface. Centrifugation makes it possible to spin down all the MT bundles onto the 
interface, leading to the formation of a dense quasi-2D MT film which exhibits local orien¬ 
tational order. Each kinesin cluster binds to multiple MTs. As each motor within the cluster 
hydrolyzes adenosine triphosphate (ATP), it moves towards the plus end of a MT and induces 
inter-filament sliding 29 . This generates extensile mechanical stresses that drive the nematic film 
away from equilibrium (Fig. lb). A biochemical regeneration system maintains a constant ATP 
concentration and powers the system for over 24 hours (see Supplementary Methods). We im¬ 
age these active nematics with both fluorescence microscopy and LC-PolScopd- 30 -. LC-PolScope 
measures the orientation of the nematic director 6(r) with pixel resolution. It also measures the 
magnitude of birefringence which reveals the effective thickness of the nematic film, or retar- 
dance A(r) (see Supplementary Figure 4 for an extended discussion). Using a 4x objective, we 
observe the dynamics of the MT film over the entire field of view, corresponding to an area of 
2.3 x 1.7 mm. 

In parallel, we have developed a tractable coarse-grained computational model. Our ap¬ 
proach employs Brownian dynamics simulations of rigid spherocylinders which, in equilibrium, 
form a nematic phased. Long-ranged hydrodynamic interactions are omitted, producing an es¬ 
sentially dry system. The length of each constituent rod increases at a constant rate, producing 
an extensile stress similar to the motor-driven extension of MT bundles (Fig. Id). Upon reach¬ 
ing a pre-set maximum length, a spherocylinder is split in half and two other rods are simul¬ 
taneously merged, thus keeping the total particle number fixed (see Supplementary Methods). 
Though inspired by the dynamics of MT bundles, this approach is not meant to quantitatively 
reproduce specific features of the experimental system, but simply to capture its microscopic 
symmetries (nematic interparticle alignment, and extensile nematic activity). 

In equilibrium, nematic defects anneal to minimize free energy, eventually producing a 
uniformly aligned state. It is not possible to prepare an analogous state in extensile active ne¬ 
matics, since uniform alignment is inherently unstable to bend deformations 5 . These distortions 
grow in amplitude and produce a fracture line, terminated at one end by a defect of charge 
+1/2, and by a —1/2 defect at the other (Fig. lc). The asymmetry of +1/2 defects causes 
motor-generated stresses to produce a net propulsive force, leading to extension of the fracture 
line. Above a critical length, the fracture line self-heals, leaving behind a pair of isolated, op¬ 
positely charged defects-! For experimental ATP concentration, +1/2 defects move at speeds 
of ~8 //m/scc. Defects are transient objects; on average, a +1/2 defect exists for 40 sec be¬ 
fore colliding with a —1/2 defect and annihilating, leaving behind a uniformly aligned nematic 
region 32 (Fig. lc). The system reaches a steady state in which the rate of defect generation is 
balanced by the rate of annihilation (Supplementary videos 1-2). Very similar patterns of defect 
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generation and annihilation are observed in the simulations, despite the fact that our computa¬ 
tional model does not include hydrodynamic interactions (Fig. le, If, Supplementary videos 
4-5). 


We have developed algorithms that identify defect positions and orientations from either 
LC-PolScope images or simulation configurations, and track their temporal dynamics over the 
entire lifetime of either a microtubule sample or simulation (see Supplemental Information). 
Defect positions and orientations are determined by measuring the winding of the local direc¬ 
tor. We define the orientation of each comet-like +1/2 defect by drawing an arrow from the 
comet’s head to its tail (see Supplemental Information). These algorithms allow us to analyze 
statistically large defect populations, providing invaluable insight into their self-organization at 
macroscopic scales. 

It has been commonly assumed that the dynamics of motile defects leads to a disordered 
chaotic state. However, quantitative analyses reveal that this is not the case. We find the that ori¬ 
entational distribution function of +1/2 defects is not flat, but exhibits two well-defined peaks, 
implying the existence of a higher-order dynamical phase (Fig. 2e). Although +1/2 defects are 
polar objects, we find that they form a nematic phase in which they are equally likely to point 
in either direction along the preferred axis (Fig. 2b). Furthermore, the orientational distribu¬ 
tion function of —1/2 defects exhibits six-fold symmetry, although the strength of this order 
is significantly less than the nematic order of +1/2 defects. Additionally, radial distribution 
functions of both ±1/2 defects reveal no long-range positional order (Supplementary Fig. 9). 

Next, we investigated how the orientational order of +1/2 defects persists in time and 
space. We find that within a single field of view (2.3 x 1.7 mm) the axis of defect order does not 
change appreciably over the entire sample lifetime (Fig. 2g). Therefore, we used a motorized 
x-y stage to repeatedly scan centimeter-sized samples every ten minutes, allowing us to measure 
long-range variations in defect ordering. In such samples we measure nearly uniform system- 
spanning orientational order (Fig. 2c, 2d). The largest active nematic sample analyzed (5cm x 
2cm) contained ~20,000 defects, demonstrating that orientational order persists at scales larger 
than 100 average defect spacings (Supplementary video 3). The defect orientational order is a 
result of spontaneously broken symmetry, and is not strongly influenced by the sample bound¬ 
aries. To demonstrate this, we have confined active nematics in a circular geometry, finding 
that defects form a single uniformly aligned domain, rather than aligning with the boundaries 
(Supplementary Fig. 1). Consistent with these observations, we additionally note that active 
nematics in rectangular channels do not strongly favor either the long or the short axis of the 
channel. 

In simulations, +1/2 defects also attain system-spanning orientational order (Fig. 3b, d). 
However, in contrast to the nematic defect ordering observed in experiments, in the computa- 
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tional system defects align with polar symmetry, which leads to their net transport along the 
preferred direction. Possible reasons for this difference are discussed below. 

Tuning filament density controls the strength of the emergent defect order and can even 
transform the system into an isotropic state. To quantify the degree of defect ordering, we mea¬ 
sure the 2D polar and nematic order parameters, P = (cos( i/j — +)) and S = (cos(2['0 — $])), 
respectively, where ^ is the orientation of a +1/2 defect and + is the mean orientation of all 
defects in a given system configuration. We find that thin nematic films (low MT concentration, 
hence low retardance) have high defect nematic order, S; increasing the film thickness (high 
MT concentration, high retardance) decreases the magnitude of S to the point where defects 
become effectively isotropic (Fig. 2a-b, 4a). A similar effect is observed in simulations when 
varying the particle density (area fraction); at the lowest densities studied, defects have rela¬ 
tively strong alignment, P. Increasing density induces a transition to an isotropic state (Figs. 
3a-b, 4b). Spatial correlation functions of these order parameters demonstrate that in all ex¬ 
perimental and computational systems with measurable defect ordering, defect correlations are 
system-spanning (Fig. 4c, 4d). Though the density of material is an easily tuneable control 
parameter, it influences many material properties, including the rate of energy dissipation, elas¬ 
tic constants, and the efficiency with which active stresses are transmitted through the material, 
all of which influence emergent properties of the system. For example, defect density in ex¬ 
periments decreases weakly with the MT film thickness, while in simulations it increases with 
rod density (Supplementary Fig. 2). Additional studies are required to relate the microscopic 
parameters to defect alignment. 

In both experiment and simulation, we find that defect-ordered systems also display ne¬ 
matic alignment of the constituent filaments (Fig. 2f, 3e). The existence of such ordering in the 
presence of large numbers of defects contrasts sharply with equilibrium systems, in which de¬ 
fects reduce or eliminate nematic order. Moreover, the ordering of both defects and constituent 
rods decreases with system density (Fig 4a-b, Supplementary Fig. 6), unlike equilibrium ly¬ 
otropic liquid crystals, in which order increases with density. These contrasts suggest that the 
orientational order in active nematics is driven by non-equilibrium dynamics, and cannot be 
accounted for purely by equilibrium-like alignment of the underlying material. 

Although +1/2 defects are preferentially created perpendicular to the local nematic field 32 } 
they do not keep this orientation but are continually reoriented by their interactions with the lo¬ 
cal environment. A +1/2 defect moving through a distorted nematic field will follow a curved 
path, always remaining perpendicular to the field in front of it, leading to complex meander¬ 
ing trajectories. Simultaneously, the passage of a +1/2 defect leaves behind a distortion in 
which the nematic field is rotated 90° from its previous orientation, creating a topological struc¬ 
ture which cannot relax except by the passage of additional defects. These distortions strongly 
affect defect motion even after the defect which formed them has moved away or been annihi- 
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lated. We therefore anticipate that a theory capable of explaining the origin of defect ordering 
will need to simultaneously account for both the defects and the underlying field. 

Since our simulation model mimics the microscopic symmetries of the experiments, it 
is notable that the two systems exhibit different emergent symmetries. In experiments, +1/2 
defects align with nematic symmetry, and the direction of defect ordering is aligned with the 
average MT direction. In simulations, defects align with polar symmetry, and the underlying 
rods are, on average, offset 90° from the defects. A number of distinctions between the two 
systems may account for these differences. First, analysis of the distributions of bend and splay 
distortions suggests that the simulation model has a much higher bend modulus (see Supplemen¬ 
tary Methods). Second, the experimental system is subjected to hard-wall boundary conditions, 
which preclude global polar ordering of motile defects. Third, the simulation model is dry, 
whereas the experiments experience hydrodynamic interactions. Finally, simulations investi¬ 
gate a true 2D system, while in experiments MTs can pass over each other. Further studies will 
be required to disentangle these effects. 

In summary, our work demonstrates that transient, short-lived motile defects can form 
higher-order dynamical phases with persistent orientational order. The existence of such phases 
in both the experimental and computational systems suggests that this is a generic feature of 
active nematics. 
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Figure 1: Overview of experimental and simulation systems, a, Microtubules (MTs) are bun¬ 
dled together by the depletion agent PEG. Kinesin clusters crosslink MTs and induce inter¬ 
filament sliding. Bundles are confined to a surfactant-stabilized oil-water interface, where they 
form a quasi-2D active nematic film, b, Fluorescence microscope image of a MT active nematic 
with defects of charge +1/2 (red) and -1/2 (blue), c, Image sequence illustrating the creation 
(top) and annihilation (bottom) of a defect pair. Scale bars 50/; m. d, Simulation microdynam¬ 
ics, consisting of hard rods which grow and split, reminiscent of the extension of MT bundles. 
L max is the length at which a rod is split, e, Snapshot of simulated active nematic with marked 
defects. Rod colors indicate their orientations, and black streamlines guide the eye over the 
coarse-grained nematic field, f, Creation (top) and annihilation (bottom) events occur analo¬ 
gously to those in experiments. Scale bars 2 L max . 
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Figure 2: Defect-ordered phase in experiments, a, Retardance map of a thick MT film in the 
regime of weak defect alignment. Red and blue markers indicate locations and orientations of 
+1/2 and —1/2 defects. Scale bar 200/im. b, Thin MT film showing strong alignment of +1/2 
defects. Scale bar 200/im. c, Orientational order in a large active nematic sample. Each red 
bar’s orientation and length indicates the mean direction and strength of defect alignment in one 
field of view. Scale bar 2mm. d, Defect alignment spans the largest samples studied (6cm x 
2cm), containing ~20,000 defects. Scale bar 10mm. e, Normalized histogram of +1/2 (red) 
and —1/2 (blue) defect orientations, P(ij)). f, MT orientation, P(9), for the sample shown in 
panels b-d. Measurements of P{6 ) from 0 to 7r are duplicated from ix to 27 t. Both P('tl') and P(9) 
show strong nematic ordering, g, The preferred defect orientation (green) and magnitude of the 
order parameter (red) averaged over a field of view persists over the entire sample lifetime. 
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Figure 3: Defect-ordered phase in simulations, a, Snapshot of a high-area-fraction simulation in 
which +1/2 defects are not aligned. Red and blue markers indicate locations and orientations of 
+1/2 and —1/2 defects. Scale bar 5L max . b, A low-area-fraction system in which defects show 
flocking-like polar alignment. Scale bar 5L max . c, A large simulation with defects aligned over 
long distances. Scale bar 5L max . d, Normalized histogram of +1/2 (red) and —1/2 (blue) defect 
orientations, P(ip). e, Rod orientations P(9) measured in the same sample. The former shows 
polar ordering, while the latter exhibits nematic ordering, and the preferred axes are offset by 
90 degrees. Measurements of P(9) from 0 to n are duplicated from 7r to 27 t. 
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Figure 4: Quantitative measurements of defect alignment, a, The defect nematic order parame¬ 
ter, S, decreases as a function of the MT film’s retardance. The blue shaded region represents 
the “noise floor” (see Supplementary Information). Inset: The nematic order parameter of the 
underlying MT filaments, S, also decreases with the MT film’s retardance. b, The polar defect 
order parameter, P, showing a transition between ordered and isotropic regimes as a func¬ 
tion of particle density. Error bars indicate a 90% confidence interval computed by bootstrap 
methods (see Supplementary Information for details), c, The two-point nematic correlation of 
defect orientation C's(r) = (cos 2 (0(r) — ”0(0))) in MT films, which shows that orientational 
order persists over long distances. Id indicates the mean inter-defect spacing. The magni¬ 
tude of ordering falls as retardance increases, d, The polar correlation of defect orientation 
Cp(r ) = (cos (0(r) — 0(0))) in simulations, a indicates the width of a single rod (see Supple¬ 
mentary Information). At short ranges, defects tend to point in opposing directions, but beyond 
the first shell of neighbors, defects are likely to be aligned in the same direction. 
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